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Abstract 

A  summary  of  the  state-of-the-art  in  nonlinear  finite  element  analysis  is 
made  by  describing  a  nonlinear  theory  and  presenting  some  case  studies.  The 
formulation  is  applicable  to  problems  of  large  displacement  and  small  strains. 

The  paper  then  focuses  on  the  general  purpose  program.  The  concept  and 
development  of  a  general  purpose  program  is  described.  A  discussion  is  then  made 
of  the  different  sizes  of  problems  which  can  be  solved  by  such  a  program.  These 
sizes  are  dependent  on  the  available  computer  core.  The  conclusion  is  made  that 
the  general  purpose  program  is  a  powerful  means  of  implementing  finite  element 
analysis  over  a  wide  spectrum  of  problems  in  structural  mechanics. 
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Introduction 

In  recent  years  we  have  seen  an  increasing  use  of  the  finite  element 
method  for  both  research  and  development.  This  paper  briefly  summarizes  the 
method  and  traces  the  reasons  for  its  widespread  use. 

The  finite  element  method  is  dependent  on  the  combination  of  two  basic 
ideas.  The  first  is  the  recognition  that  problems  in  continuum  mechanics  may  be 
solved  by  complete  satisfaction  of  only  one  of  the  two  requirements  of  equilibrium 
and  compatibility  if  the  other  condition  is  also  satisfied  in  an  integral  sense. 
This  approximate  solution  of  the  remaining  condition  in  the  integral  sense  is 
brought  about  by  the  use  of  the  principle  of  virtual  work  and  theorems  resulting 
from  it.  The  second  idea  is  that  the  function  for  a  whole  domain  may  be  better 
approximated  by  local  functions  assumed  within  subdomains  which  also  maintain  con¬ 
tinuity  of  the  functions  across  the  subdomains.  The  undetermined  parameters  for 
the  assumed  local  subdomain  functions  are  then  related  to  physical  quantities  of 
displacement  [1]  or  force  [2]  degrees  of  freedom  at  points  or  nodes  on  the  bound¬ 
aries  of  the  subdomain.  This  then  allows  the  definition  of  equations  which  define 
either  stiffness  or  flexibility  matrices  for  the  subdomain  (or  discrete  element). 
The  combination  of  the  relaxation  of  the  requirement  of  either  equilibrium  or  com¬ 
patibility  and  the  localized  functions  whose  unknowns  are  represented  by  quan¬ 
tities  at  nodes  results  in  a  considerable  easing  of  the  problems  of  geometry  and 
boundary  conditions.  Different  subdomains  may  be  modeled  with  different  functions 
and  these  may  be  used  simultaneously  for  an  analysis. 

The  finite  element  theory  is  usually  cast  in  matrix  theory  since  this 
allows  the  large  background  of  matrix  theory  to  be  exploited.  Its  development 
occurred  at  the  same  time  as  the  order  of  magnitude  increase  in  computer  speeds 
and  core  size.  This  happy  confluence  of  all  the  factors  discussed  above  has 
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given  rise  to  the  widespread  development  and  use  of  the  finite  element  method. 
Initially,  its  implementation  took  place  in  the  form  of  specialized  programs 
written  for  specific  purposes.  Then  as  the  method  developed  it  became  obvious 
that  a  more  general  approach  could  be  adopted  in  which  the  common  tasks  for  every 
finite  element  could  be  programmed  once  and  for  all.  This  has  resulted  in  the 
development  of  the  general  purpose  finite  element  programs.  These  efforts  have 
the  same  overall  strategy  as  the  SPADE  projects  adopted  for  partial  differential 
equations.  At  present,  the  continued  development  of  the  general  purpose  programs 
appears  to  be  the  best  means  of  implementing  finite  element  theory.  Yet  notwith¬ 
standing  this ,  little  has  appeared  in  the  literature  which  specifically  concerns 
itself  with  the  features  and  underlying  philosophy  of  the  general  purpose  program. 
It  is  the  purpose  of  the  current  paper  to  discuss  the  development  of  a  general 
purpose  program. 

Review  of  Literature 

The  present  paper  will  trace  developments  from  the  original  paper  by 
Turner  et  al.  [1],  and  attention  will  be  confined  solely  to  the  direct  stiffness 
method  of  finite  element  analysis. 

Initially,  work  was  concerned  with  developing  elements  [3-8]  with  compat¬ 
ible  displacements  at  the  boundaries.  This  phase  can  now  be  said  to  be  complete, 
and  elements  exist  to  cover  any  two-  or  three-dimensional  solids  (including  shell 
structures).  We  may  classify  elements  by  the  type  of  displacement  modes  assumed 
and  with  this  classification  three  types  of  elements  can  be  recognized.  In  its 
two-dimensional  form  these  three  may  be  referred  to  as  the  triangular,  the  orthog¬ 
onal  and  the  piecewise  patching  type  of  element. 

In  the  triangular  type  of  element,  the  displacement  modes  are  assumed  to 
take  the  form  of  complete  polynomials  [1-4].  In  the  orthogonal  type  of  element. 
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the  displacement  modes  are  assumed  to  take  the  form  of  either  Lagrange  Polynomials 
[3,5]  or  Hermitian  Polynomials  [6].  The  Lagrange  Polynomials  [5]  are  also  used 
to  extend  the  element  formulation  to  quadrilaterals  and  curve  shaped  elements  by 
isoparametric  techniques.  In  the  patching  type  of  element  [7,8],  usually  used  for 
shells  and  referred  to  as  the  DeVeubeke  element,  the  displacement  modes  are 
assumed  to  be  made  up  of  a  compatible  patching  of  complete  polynomials.  Three- 
dimensional  equivalents  also  exist  for  the  first  two  types  of  elements. 

At  the  same  time  work  [9,10]  was  reported  which  established  a  framework 
by  which  the  finite  element  method  could  be  related  to  the  methods  used  in  contin¬ 
uum  mechanics.  The  finite  element  method  is  now  recognized  as  a  special  case  of 
the  Rayleigh-Ritz  method  where  generalized  modes  are  assumed  over  subdomains  where 
the  generalized  modes  give  rise  to  inter-element  compatibility  of  displacements. 

With  the  establishment  of  the  method,  attention  was  turned  to  extensions 
for  nonlinear  analysis.  In  the  area  of  material  nonlinearity,  two  methods  were 
developed  for  elastic-plastic  analysis.  The  method  of  initial  strains  is  based 
on  the  idea  of  modifying  the  equations  of  equilibrium  so  that  the  elastic  equa¬ 
tions  can  be  used  throughout  on  the  left-hand  side  of  the  equations.  Modifica¬ 
tions  are  introduced  on  the  right-hand  side  of  the  equation  to  compensate  for  the 
fact  that  the  plastic  strains  do  not  cause  any  change  in  the  stresses.  On  the 
other  hand,  the  tangent  modulus  method  is  based  on  the  linearity  of  the  incremen¬ 
tal  laws  of  plasticity  and  approaches  the  problem  in  a  piecewise  linear  fashion. 
The  load  is  applied  in  increments ,  and  at  each  increment  a  new  set  of  coeffi¬ 
cients  is  obtained  for  the  equilibrium  equations.  The  matrix  equations  for  the 
finite  element  analysis  using  the  method  of  initial  strains  were  developed  by 
Padlog  et  al.  [11],  Argyris  et  al.  [12]  and  Jensen  et  al.  [13].  A  recent  paper 
by  Witmer  [14]  summarizes  the  latest  application  of  the  method.  The  equations 
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for  the  tangent  modulus  method  were  developed  by  Pope  [15],  Swedlow  and  Yang  [16] 
and  Marcal  and  King  [17].  The  two  methods  were  compared  by  Marcal  [18]  and  a 
close  similarity  was  found  between  them. 

Progress  has  also  been  made  in  the  area  of  geometric  nonlinearity.  Large 
displacement  analysis  by  the  finite  element  method  was  first  proposed  by  Turner 
et  al.  [19].  Initial  stress  stiffness  matrices  were  developed  to  account  for  the 
effect  of  initial  stress  in  truss  and  plane  stress  assemblies.  Subsequent  work 
on  the  derivation  of  the  initial  stress  matrices  for  other  elements  were  reported 
by  Argyris  et  al.  [20],  Gallagher  et  al.  [21]  and  Kapur  and  Hartz  [22].  Martin 
[23]  placed  the  derivation  of  the  initial  stress  matrix  on  a  firm  foundation  by 
using  a  potential  energy  formulation  together  with  the  nonlinear  strain  displace¬ 
ment  relation  (for  Green's  strain).  The  above  papers  were  concerned  with  forming 
matrices  which  account  for  geometric  changes  in  the  solid  during  an  increment  of 
load.  These  matrices  were  then  either  used  in  a  piecewise  linear  manner  or  used 
in  an  eigenvalue  analysis  of  the  Euler  type.  Recent  papers  [24,25]  have  drawn 
attention  to  the  fact  that  certain  important  terms  were  neglected  in  finite  ele¬ 
ment  large  displacement  analysis.  These  neglected  terms  result  in  what  was  called 
the  initial  displacement  matrix  and  is  a  result  of  the  coupling  between  the  quad¬ 
ratic  and  the  linear  terms  in  the  strain  displacement  expressions. 

Other  workers  solved  the  nonlinear  equations  of  the  finite  element  method 
directly.  Bogner  et  al.  [26]  performed  a  direct  minimization  of  the  potential 
energy  without  explicitly  forming  the  matrix  stiffness  equations.  The  large  dis¬ 
placement  behavior  was  followed  into  the  post-buckling  region.  Mallett  and 
Berke  [27]  applied  this  method  to  frameworks  and  Bogner  et  al.  to  plates  and 
shells  [28].  Oden  [29]  and  Oden  and  Kubitza  [30]  developed  nonlinear  stiffness 
relations  for  the  nonlinear  elasticity  problem.  The  equations  were  solved  by  a 
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Newton-Raphson  method.  This  series  of  works  is  perhaps  best  placed  in  a  separate 
category.  It  concerns  itself  with  large  strain,  large  displacement  analysis.  The 
other  papers  reviewed  previously  have  all  implicitly  assumed  a  large  displacement 
small  strain  theory.  In  addition,  the  constitutive  equations  used  there  are  in 
terms  of  an  energy  potential  which  is  appropriate  for  a  rubber-like  material.  A 
similar  formulation  with  appropriate  assumptions  of  constitutive  behavior  gives 
rise  to  equations  for  large  strain  large  displacement  analysis  of  metal  structures 
[31].  A  recent  survey  by  Oden  [32]  brings  out  the  very  general  nature  of  the 
finite  element  formulation. 

Hence,  we  have  seen  that  progress  has  been  made  in  both  nonlinear  material 
and  geometric  behavior.  The  two  nonlinear  formulations  do  not  depend  on  each 
other  so  that  they  may  be  profitably  combined.  In  the  present  paper  we  shall 
restrict  our  attention  to  a  small  strain  large  displacement  theory  appropriate  to 
shells  and  other  solid  metal  structures. 

In  the  area  of  elastic  analysis  by  the  finite  element  method,  general  pur¬ 
pose  computer  programs  exist  which  are  written  with  a  view  to  covering  the  whole 
area  of  stress  analysis.  These  general  purpose  programs  exploit  the  generality 
of  the  matrix  formulation  of  the  finite  element  method.  The  programs  have  a 
library  of  elements  which  can  be  used  for  the  modeling  of  most  structures  in 
service.  Melosh  et  al.  [33]  have  summarized  the  more  recent  general  purpose  pro¬ 
grams  . 

Technical  Considerations 

The  theory  outlined  here  has  been  developed  previously  in  [34],  It  is 
included  here  for  completeness.  The  displacement  method  of  finite  element  analysis 
will  be  used  throughout.  The  structure  to  be  analyzed  is  divided  into  a  number 
of  elements.  The  behavior  of  each  element  is  lumped  into  a  number  of  nodal  point 


6 


displacements.  Conforming  displacement  modes  or  simply  conforming  elements  are 
modes  which  maintain  displacement  compatibility  between  adjacent  elements  at  the 
element  boundaries.  The  principle  of  virtual  work  is  then  used  to  effect  the 
lumping  of  the  equivalent  forces.  The  principle  of  virtual  work  is  of  course 
applicable  to  large  displacement  as  well  as  nonlinear  material  behavior. 

A  brief  outline  of  the  method  is  now  given .  A  displacement  mode  is  first 
chosen  for  the  type  of  element  being  used, 

n 

u  =  y  a.f.(x)  =  [f(x)]{a}  (1) 

i=l  1  1 

where  u  is  the  displacement  at  position  x 

x  is  used  to  represent  the  coordinates  of  the  element 
a^  are  the  generalized  displacements  (also  written  {a}) 
n  is  the  number  of  terms  in  the  summation. 

By  substituting  for  x  at  the  nodes  obtain 

{a}  =  Ca]{u}  (2) 

where  {u}  is  the  displacement  at  the  nodal  points  (note  that  there  is  a  dis¬ 
tinction  between  the  bracketed  and  unbracketed  u  ). 

[a]  is  the  nodal  point  to  generalized  displacement  transformation  matrix. 
By  the  assumption  of  small  strains  we  also  have 

A{a}  =  [a]A{u}  (3) 

where  the  prefix  A  denotes  an  increment  of  the  quantity  immediately  following 
it. 

We  now  define  the  so-called  differential  operator  [B]  which  transforms 
an  increment  of  generalized  displacement  to  an  increment  of  strain. 


(4) 
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Me)  =  [B]A{a) 

The  differential  operator  [B]  is  a  function  of  position  x  and  of 
current  displacement  u  .  It  is  defined  by  writing  the  nonlinear  strain  dis¬ 
placement  equations  (Green's  strain)  in  incremental  form. 

The  strain  increment  can  be  written  as  a  stress  increment  for  an  elastic- 
plastic  material  in  the  manner  of  Marcal  and  King  for  solids  [17]  and  Marcal  T35] 
for  plates  and  shells. 

A {a}  =  [D]A{e>  (5) 

The  stress  increments  A{o)  are  accumulated  at  representative  points  within  each 
element.  Each  representative  point  is  given  a  set  of  reference  axes  which  deform 
with  the  element  and  so  take  the  same  direction  as  that  defined  by  an  increment 
of  Green's  strain.  Thus  stresses  and  strain  increments  are  automatically  aligned 
and  the  nonlinear  equations  of  equilibrium  can  be  set  up  with  ease. 

We  now  use  the  principle  of  virtual  work  to  define  equivalent  forces  {P} 
at  the  nodes  for  a  virtual  displacement  6 [uj  . 

A[uJ{P}  =  [  5[eJ{o}dV  =  [  A  |_u|  [ct]T[B]T{a}dV  (6) 

Jy  Jy 

where  [_  J.  denotes  a  row  vector,  and  integration  is  performed  over  the  volume  V 
Cancelling  the  non-zero  virtual  displacements  from  both  sides  and  writing 
equation  (6)  in  incremental  form  with  the  aid  of  equation  (5),  obtain 

A{P}  =  [  [a]TA[B]T{a)dV  +  f  [a]T[B]T[D][B][a]dVA{u}  (7) 

Jy  Jy 

The  last  term  on  the  right  of  equation  (7)  can  be  divided  into  a  matrix  which  is 
dependent  on  the  current  displacement  and  one  which  is  not.  With  some  rearrange¬ 
ment  ,  we  obtain  the  element  stiffness  matrices  <*,. 

j  rv  t  t.xt 

Fl'.j ‘/TT7."  fir.  , 

fi’A’ISAP-jx 
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A  { P }  =  ([kU)]  +  [k(2)]  +  [k(0)])  A{u} 


(8) 


where  [k*'1^]  is  the  initial  stress  matrix  and  is  obtained  from  the  first  term 
on  the  right  of  equation  (7). 

(2) 

[k  ]  is  the  initial  displacement  matrix  of  Marcal  [25]. 

[k^^]  is  the  small  displacement  stiffness  matrix. 

The  element  stiffness  matrices  and  the  nodal  equivalent  forces  are  then 
summed  in  the  usual  direct  stiffness  manner  to  obtain  master  stiffness  equations 
represented  by  equations  in  capitals 

,(l)n  .  r./(  2 )  -,  .  ri,(0). 


A{P}  =  ([k'1']  +  [K^;]  +  [K^U  j) A{u} 


(9) 


General  Purpose  Programs 

The  stiffness  equations  developed  are  quite  general  and  are  not  restricted 
to  a  particular  type  of  element.  It  is  therefore  possible  to  write  a  general  pur¬ 
pose  program  which  implements  the  theory.  This  program  will  then  serve  as  a  basic 
and  common  program  to  which  subroutines  may  be  added  to  account  for  specific  char¬ 
acteristics  belonging  to  the  particular  type  of  element  (or  element  combinations) 
being  used.  Two  general  approaches  to  programming  language  have  been  adopted. 

The  one  most  favored  by  the  developers  of  ASKA,  FORMAT,  MAGIC,  NASTRAN,  SAMIS  and 
STRUDL  is  to  make  use  of  some  type  of  matrix  interpretive  language.  Here  the 
intention  is  to  develop  machine  independent  concepts  and  also  lay  the  foundation 
for  easy  implementation  of  further  theoretical  developments.  However,  most  of 
these  programs  have  been  developed  under  the  influence  of  particular  computers  and 
programming  languages,  so  that  the  aims  of  complete  machine  independence  and  free¬ 
dom  from  bookkeeping  requirements  have  not  been  fully  achieved.  On  the  other 
hand,  there  have  been  other  programs  (ELAS,  MARC2)  which  were  developed  with 
FORTRAN  as  the  programming  language  and  making  use  of  matrix  support  packages. 
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The  one  common  feature  to  both  approaches  is  the  attempt  to  implement  the 
matrix  manipulation  required  by  the  theory  in  as  general  a  form  as  possible.  We 
see  from  (7)  and  (8)  that  the  matrix  operations  required  to  form  the  stiffness 
matrices  do  not  change.  Similarly,  the  assembly  of  the  element  stiffness  to  form 
the  master  stiffness  matrix  does  not  change  with  different  elements. 

In  order  to  focus  on  the  advantages  of  general  purpose  programs,  we  shall 
now  focus  on  the  program  MARC2  developed  at  Brown  University.  Most  of  the  features 
found  in  this  program  can  readily  be  included  in  other  programs  so  that  the  points 
to  be  discussed  can  be  thought  to  apply  equally  to  all  programs  in  general.  This 
program  was  developed  with  the  intention  of  carrying  out  the  common  matrix  opera¬ 
tions  required  to  solve  finite  element  problems  with  nonlinear  material  and/or 
geometric  behavior.  Because  it  was  meant  to  be  used  in  a  research  environment, 
it  was  organized  with  a  view  to  minimize  the  coding  required  to  implement  new 
elements.  This  is  made  easier  by  the  use  of  numerical  integration  to  form  the 
element  stiffness  matrices.  Only  four  user  subroutines  are  required  to  form  the 
[<*][B]  quantities  and  specify  the  weighting  functions  required  to  perform  the 
numerical  integration.  The  general  purpose  program  carries  out  the  rest  of  the 
calculations  based  on  input  data.  In  particular,  a  subroutine  has  been  developed 
to  implement  the  incremental  Prandtl-Reuss  relations.  Another  subroutine  inte¬ 
grates  these  relations  through  the  thickness  for  a  plate  or  shell  when  required. 
Various  subroutines  enable  the  assembled  matrix  equations  to  be  solved  by  either 
the  direct  or  iterative  approach,  as  well  as  giving  the  option  of  an  in-core 
assembly  and  out-of-core  solution.  This  program  simplifies  combined  elastic- 
plastic  or  creep  and  large  displacement  analysis  by  reducing  the  amount  of  addi¬ 
tional  programming  required  from  a  user.  The  nonlinear  problem  is  converted  to 
a  series  of  piecewise  linear  problems.  There  is  now  an  increase  of  an  order  of 
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magnitude  in  computing  time  required  compared  to  a  linear  elastic  solution  since 
it  takes  about  ten  steps  to  trace  the  load  history  of  a  structure  to  its  buckling 
or  limit  load. 

Figure  1  shows  the  flow  chart  for  the  general  purpose  program  MARC2. 

The  procedures  depicted  in  the  main  flow  are  the  control,  assembly,  application 
of  boundary  conditions  and  the  solution  of  the  master  stiffness  equation.  Two 
subroutines  interface  with  the  user  subroutines  and  form  the  element  stiffnesses 
and  the  initial  stress  and  strain  vectors  respectively.  Their  purpose  is  to 
organize  and  perform  the  numerical  integration  to  obtain  the  required  quantities. 
In  turn,  these  subroutines  draw  on  the  subroutines  which  form  the  linear  incre¬ 
mental  strain  to  stress  transformation  matrix  [D]  referred  to  above. 

It  is  of  interest  to  note  here  the  various  types  of  problems  that  can  be 
handled  by  the  program.  It  is  noted  that  these  can  be  performed  with  any  com¬ 
binations  of  elements  and  any  combination  of  the  following  classes  of  problems : 

1.  Elastic 

2.  Elastic-plastic 

3 .  Creep 

4 .  Thermal  strains 

5 .  Large  displacement 

6 .  Large  strains 

7.  Buckling  (eigenvalue  analysis  at  any  load  level) 

We  see  immediately  the  advantages  of  using  a  general  purpose  program.  Any 
feature  implemented  in  the  program  can  be  combined  with  all  previous  developments. 
As  illustrations  of  this  we  give  examples  of  two  recent  additions  to  the  general 
purpose  program.  The  first  was  the  implementation  of  a  buckling  analysis.  Once 
this  feature  was  checked  it  meant  that  it  was  possible  to  make  use  of  all 
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previously  developed  elements  with  large  displacement  capabilities  and  perform 
buckling  analysis  of  beams,  plates  and  shells.  Conversely,  as  an  example  of  using 
the  common  features  in  the  program,  a  new  arbitrary,  doubly  curved  shell  element 
[36]  was  recently  developed.  It  was  then  possible  to  use  the  element  in  the  solu¬ 
tion  of  all  the  seven  classes  of  problems  outlined  above.  This  ability  to  pre¬ 
serve  and  exploit  all  previous  developments  is  the  main  reason  behind  the  impetus 
towards  development  of  general  purpose  programs  in  recent  years.  This  generaliza¬ 
tion  is  in  accord  with  the  development  and  use  of  computer  programs  in  other  areas. 
One  other  advantage  of  the  general  purpose  program  is  that  the  main  flow  of  this 
program  will  be  used  frequently  and  confidence  in  such  a  program  will  be  more 
readily  established. 

There  are  also  some  drawbacks  in  such  a  general  approach  which  are  perhaps 
not  so  evident.  First  of  all  there  is  its  slower  running  time  because  of  the  many 
conditional  statements  in  the  program.  This  slower  speed  is  particularly  notice¬ 
able  in  the  larger  computer  systems  where  parallel  computing  devices  are  employed. 
Such  a  program  also  tends  to  become  large,  particularly  if  there  is  a  large  team 
working  on  it,  and  the  limits  of  computer  storage  are  quickly  reached.  Another 
problem  to  be  overcome  is  that  of  verification  documentation  and  dissemination  of 
such  a  program.  Because  MARC2  was  intended  to  be  used  in  a  research  environment, 
the  coding  has  been  kept  to  a  minimum.  Even  so,  it  has  grown  to  about  6000  FORTRAN 
statements  and  already  makes  severe  demands  on  new  users.  It  is  interesting  to 
note  here  that  it  takes  a  new  user  about  a  month  and  a  half  to  learn  the  program 
and  begin  to  contribute  to  its  development  by  modifying  it.  One  major  disadvantage 
in  developing  such  a  program  is  the  difficulty  in  keeping  changes  made  by  one 
worker  from  interfering  with  the  programming  work  of  others. 
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Mote  on  the  Cost  of  Computing 

In  this  note  we  shall  examine  the  relationship  between  the  size  of  a 
problem  and  the  cost  of  computing.  The  actual  cost  of  computing  depends  on  the 
system  configuration  and  it  is  possible  to  obtain  differences  in  costs  of  up  to 
factors  of  2  by  simply  choosing  different  machines  of  the  same  nominal  speed,  as 
well  as  by  choosing  the  same  make  of  computer  but  using  it  at  different  installa¬ 
tions.  The  important  point  to  be  recognized  is  that  core  is  now  required  on  the 
faster  machines  in  such  large  sizes  that  its  cost  is  as  much  as  that  of  the 
central  processor  (C.P.U.).  Thus  realistic  accounting  procedures  recognize  this 
and,  since  particular  system  configurations  are  designed  with  certain  functions 
in  mind,  its  use  for  other  purposes  may  cause  a  certain  penalty.  We  shall  use 
here  the  total  system  time  as  a  measure  of  the  computing  required  to  solve  each 
problem  and  not  merely  the  C.P.U.  time  used.  The  size  of  a  problem  is  dependent 
on  the  system  configuration  and  will  differ  from  one  machine  to  the  next.  Thus 
the  following  discussion  is  an  attempt  to  measure  the  relative  cost  of  solving 
relatively  large  to  small  size  problems  on  a  particular  computer.  The  ability  to 
handle  large  problems  is  dependent  on  the  size  of  core  available  to  the  user  and, 
at  the  same  time,  it  is  also  dependent  on  the  core  left  for  simultaneous  use  by 
other  users.  Because  the  core  requirements  for  a  finite  element  problem  are 
determined  by  its  master  stiffness  equation,  we  shall  measure  the  size  of  the 
problem  by  the  product  of  the  number  of  degrees  of  freedom  with  the  half-band- 
width.  With  this  definition  the  size  of  a  problem  can  be  divided  into  three 
distinct  categories  which  are  determined  by  the  peripheral  storage  required  to 
solve  the  master  stiffness  equations.  These  solutions  fall  into  the  following 
categories,  viz.,  the  in-core  assembly  and  solution,  in-core  assembly  and 
out-of-core  solution,  and  the  out-of-core  assembly  and  solution.  In  the  first 
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category  the  complete  solution  can  be  effected  in-core.  In  the  second,  the  core 
requirements  are  such  that  assembly  of  the  matrix  can  be  performed  in  core  by 
packing  the  matrix  while  the  solution  is  carried  out  with  the  aid  of  magnetic 
discs  or  tapes.  This  course  of  action  usually  doubles  the  size  that  can  be 
handled  by  an  in-core  solution.  Finally,  the  third  category  is  one  in  which  the 
problem  is  so  large  that  both  assembly  and  solution  must  be  performed  out  of  core. 
Figure  2  shows  a  plot  of  total  system  time  used  against  the  size  of  problem.  The 
limits  of  the  in-core  solution  and  the  partial  in-core  solution  are  shown. 

Because  of  the  relative  speeds  of  C.P.U.  to  I./O.  operations  the  total  system  time 
is  shown  increasing  with  increasing  computer  speeds.  No  attempt  is  made  to  show 
relative  computing  costs.  Figure  2  also  shows  a  curve  for  computing  on  a  time- 
shared  machine  which  is  able  to  simulate  practically  unlimited  core  size  for  the 
user,  the  so-called  virtual  machines.  It  appears  that  operation  on  such  a  com¬ 
puter  does  not  penalize  a  user  unduly  for  working  out  of  core  since  this  is  the 
normal  mode  of  operation  for  which  the  system  is  planned.  A  distinct  advantage 
in  cost  can  be  gained  by  solving  the  larger  problems  on  such  a  machine.  The 
writer's  experience  does  bear  this  out. 

The  demands  on  the  computer  is  not  the  only  cost  involved  in  a  finite  ele¬ 
ment  analysis  since  the  effort  required  for  coding  and  the  preparation  of  data 
must  also  be  taken  into  account.  The  question  of  overall  economics  deserves 
further  attention.  It  would  be  interesting,  too,  to  combine  this  with  further 
study  of  future  hardware  and  software  development  of  computer  systems. 

Case  Studies 

In  this  section  we  present  a  series  of  case  studies  which  illustrate 
various  facets  of  the  current  state-of-the-art  of  finite  element  analysis.  It 
is  hoped  that  these  studies  taken  together  will  give  an  overall  picture  of  the 
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progress  made  in  this  area.  The  writer  has  mainly  drawn  from  results  obtained  in 
conjunction  with  his  colleagues.  Other  choices  could  also  have  been  made,  but 
the  present  ones  simply  reflect  greater  familiarity  with  the  results,  as  well  as 
ready  access  to  it. 

1.  Substructural  Analysis  of  the  747  Aircraft  Wing-Body  Intersection  [37] 

This  example  of  an  elastic  analysis  is  included  to  show  the  large-scale 
problems  that  are  handled  by  the  finite  element  method.  The  method  of  substruc¬ 
tures  or  matrix  partitioning  is  found  to  be  the  best  way  of  reducing  the  problem 
to  manageable  proportions  in  both  the  data  handling  and  the  equation  solving 
aspects  of  the  problem.  The  substructures  are  shown  in  schematic  form  in  Fig.  3. 
These  were  idealized  by  a  combination  of  rod,  beam,  shear  and  constant  strain 
elements.  The  whole  problem  resulted  in  13,870  degrees  of  freedom.  The  sub¬ 
structuring  reduced  the  largest  band-width  that  had  to  be  handled  at  any  one  stage. 
In  connecting  the  substructures  there  were  a  total  of  709  degrees  of  freedom  that 
interacted  at  the  interface.  The  effort  that  is  involved  in  performing  the  anal¬ 
ysis  of  the  problem  is  described  in  [37].  The  problem  is  restricted  to  linear 
elastic  behavior;  however,  with  the  current  rate  of  progress  in  the  area,  it  is 
not  difficult  to  envisage  the  same  problem  being  solved  with  nonlinear  material 
and  geometric  behavior. 

It  is  of  interest  to  note  that  about  a  hundred  man-months  of  effort 
stretching  over  seven  months  was  required.  Much  of  the  model  idealization  and 
work  on  the  substructures  proceeded  in  parallel.  Twenty-eight  hours  of  CDC  6600 
C.P.U.  time  and  one  hundred  and  twenty  hours  of  residency  time  was  required  for 
an  error- free  pass  through  the  system. 
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2.  Elastic-Plastic  Analysis  of  Tensile  Specimen  with  Semi-Elliptic  Crack  [38] 

The  next  example  is  one  of  the  analysis  of  a  semi-elliptic  crack  in  a 
tensile  specimen  using  a  combination  of  648  three-dimensional  cubic  and  isopara¬ 
metric  elements  [38]  and  900  nodal  points  or  2,700  degrees  of  freedom.  Eight 
layers  of  elements  were  used  to  model  the  problem.  Figure  4  shows  the  bottom 
layer  of  elements  and  the  layout  adopted  to  represent  the  semi-elliptic  crack. 

The  solution  was  carried  out  by  first  obtaining  the  load  to  cause  the  most  highly 
stressed  element  to  yield.  Four  increments,  each  equal  to  0.07  of  the  elastic 
load  were  added  to  study  the  elastic-plastic  behavior.  The  progress  of  plastic 
yielding  in  the  second  and  fourth  load  increment  is  shown  in  Fig.  5.  The  elastic 
solution  took  45  minutes  of  C.P.U.  time  on  the  IBM  360-91,  and  subsequent  elastic- 
plastic  increments  took  about  15  minutes  per  increment.  A  method  of  successive 
over- relaxation  was  used. 

3.  Analysis  of  Shell-Nozzle  Junction  with  Combined  Shell  and  Triangular  Ring 
Elements  [39] 

This  example  is  included  to  show  a  combination  of  shell  and  solid  elements 
by  the  method  of  linear  constraints  [39].  A  mild  steel  shell  nozzle  junction 
under  pressure  was  studied  experimentally  by  Dinno  and  Gill  [40].  This  same  prob¬ 
lem  was  analyzed  using  the  mesh  in  Fig.  6.  Triangular  ring  elements  are  used  in 
and  around  the  weld  section,  and  shell  elements  are  used  throughout  the  main  body 
of  the  shell  and  nozzle.  Comparison  of  the  finite  element  results  with  experi¬ 
mental  data  is  shown  in  Fig.  4.  The  actual  differences  between  the  peak  stresses 
can  be  seen  in  Table  1.  The  hybrid  finite  element  results  show  considerable 
improvement  over  a  previous  modified  shell  theory  approach  using  a  band  of  pres¬ 
sure  for  the  junction  [41].  That  theory  was  itself  a  large  improvement  over  the 
simple  shell  theory. 
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LIMIT  OF  PROPORTIONALITY 
(lb/ in2 

Experimental  [40] 

800 

Simple  shell  theory  [41] 

340 

Band  theory  [41] 

630 

Hybrid  analysis 

793 

TABLE  1.  First  Yield  of  Shell-Nozzle  Junction 
with  Internal  Pressure 

4.  Imperfect  Hemisphere  under  External  Pressure  [42] 

This  example  shows  the  combined  effect  of  nonlinear  geometric  and  material 
behavior.  Both  of  these  act  together  to  drastically  weaken  the  load  resistance 
of  the  structure.  The  oblate  spherical  shell  of  Fig.  8  was  analyzed  by  Bushnell 
[43]  in  the  nonlinear  elastic  region.  It  was  there  observed  that  high  "elastic" 
stresses  were  observed  at  the  crown  of  the  shell  prior  to  collapse.  This  shell 
was  analyzed  with  a  dimensionless  yield  stress  of  0.00666  E  together  with  a  linear 
work -hardening  curve  with  a  slope  of  0.05  E.  This  corresponds  roughly  to  the 
stress  strain  curve  of  an  Aluminum  Alloy. 

Figure  9  gives  a  comparison  between  the  buckling  pressure  of  the  elastic 
shell  of  [43]  and  the  present  elastic-plastic  results.  The  results  are  plotted 
in  terms  of  the  parameters  used  in  [43].  The  classical  buckling  load  pc  is 
defined  by 

p  =  1.21(2H/R)2E 
c 

where  E  is  the  Young's  Modulus 

R  is  the  radius  of  the  sphere 
and  H  is  the  half- thickness . 
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The  geometric  parameter  A  is  defined  by 

n  i/2 

A  =  /12(l-v2)  (— )  F— 

imp. 

where  R.  =  mean  radius  of  the  oblate  Dortion  of  the  sphere,  v  =  Poisson's 
imp . 

Ratio. 

Plastic  yielding  has  a  considerable  effect  on  the  behavior  of  the  oblate 
shells  under  external  pressure.  This  effect  increases  with  the  thickness  to 
radius  ratio.  It  is  noted  that  the  failures  at  the  higher  thickness  to  radius 
ratios  (A  £  1.5)  are  due  to  membrane  yield. 

Reasonable  agreement  is  also  obtained  with  the  elastic  results  of  Bushnell 
[43]  for  the  thinner  shells  which  do  not  yield  before  buckling. 

5.  Infinite  Incompressible  Log  [44] 

The  infinite  log  under  symmetric  line  loading  is  shown  in  Fig.  10.  This 
problem  was  solved  by  Oden  [44]  using  the  generalized  Newton-Raphson  method.  The 
problem  was  reduced  to  22  simultaneous  equations  by  taking  advantage  of  symmetry. 
For  illustration  purposes,  the  line  loading  P  was  taken  to  be  200  lb/in  and 
for  the  Mooney  constants  =  43.75  lb/in2  and  C ^  ~  6.25  lb/in2  were  used. 
Results  in  the  form  of  the  deformed  profile  are  indicated  in  Fig.  11. 

Conclusions  and  Future  Work 

In  this  paper  we  have  examined  the  formulation  and  the  implementation  of 
a  theory  for  nonlinear  finite  element  analysis.  The  general  purpose  program  was 
shown  to  be  a  versatile  and  flexible  method  of  implementing  the  basic  theory. 

It  was  found  possible  to  classify  three  basic  sizes  of  problems  which  were 
dependent  on  the  ability  of  the  computer  to  either  assemble  or  solve  the  master 


stiffness  matrix  in  core. 
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Case  studies  were  given  to  illustrate  representative  applications  of  the 
theory.  If,  as  is  argued  here,  the  general  purpose  program  is  the  key  to  wide¬ 
spread  applications  of  the  theory,  then  much  more  remains  to  be  learnt  about  its 
development  and  organization.  Little  is  known  about  the  best  way  to  match  pro¬ 
grams  with  particular  computer  system  configurations.  Even  less  is  known  about 
the  impact  of  future  hardware  developments.  Finally,  rigorous  procedures  have 
yet  to  be  developed  for  verification  of  these  programs. 
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FIG.  3  SCHEMATIC  OF  SUBSTRUCTURES  FOR  BOEING  747 


ELEMENTS  IDENTIFICATION  IN  THE  CRACK  PLANE 


1)  PLASTIC  ZONE  ABCDEF  IS  AFTER  TWO  INCREMENTS 

2)  PLASTIC  ZONE  A  B' B"C'C" D'D " F'F "  E'E  IS  AFTER 
4th  INCREMENT 


FIG. 5  PROGRESS  IN  PLASTIC  Yl  ELDING  TENSILE  PLATE 
WITH  SEMI-ELLIPTIC  CRACK 
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FIG. 7  NOZZLE  DISPLACEMENT 


R  imp  /  R  =  1.15 
a  =  10° 


FIG.  8  EXTERNALLY  PRESSURIZED  IMPERFECT 
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